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Abstract 

This technical report documents the theoretical, computational, and practical as- 
pects of the one-dimensional Navier-Stokes finite element flow model. The docu- 
ment is particularly useful to those who are interested in implementing, validating 
and utilizing this relatively-simple and widely-used model. 

Keywords: one-dimensional flow; Navier-Stokes; Newtonian fluid; flnite ele- 
ment; elastic vessel; interconnected network; blood flow in large vessels; branching 
flow; time-independent flow; time-dependent flow. 
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1 Introduction 

The one-dimensional (ID) Navier-Stokes flow model in its analytic formulation and 
numeric implementation is widely used for calculating and simulating the flow of 
Newtonian fluids in large vessels and in interconnected networks of such vessels [1- 
5]. In particular, the model is commonly used by bioengineers to analyze blood flow 
in the arteries and veins. This model can be easily implemented using a numeric 
meshing technique, such as finite element method, to provide a computational 
framework for flow simulation in large tubes. The model can also be coupled with 
a pressure-area constitutive relation and hence be extended to elastic vessels and 
networks of elastic vessels. Despite its simplicity, the model is reliable within its 
validity domain and hence it can provide an attractive alternative to the more 
complex and costly multi-dimensional flow models in some cases of flow in regular 
geometries with obvious symmetry. 

The roots of the ID flow model may be traced back to the days of Euler who ap- 
parently laid down its mathematical foundations. In the recent years, the ID model 
became increasingly popular, especially in the hemodynamics modeling. This is 
manifested by the fact that several researchers [2-4, 6-18] have used this model 
recently in their modeling and simulation work. 

The 'ID' label attached to this model stems from the fact that the 6 and r 
dependencies of a cylindrically-coordinated vessel are neglected due to the axi- 
symmetric flow assumption and the simplified consideration of the flow profile 
within a lumped parameter called the momentum correction factor. Therefore, 
the only dependency that is explicitly accounted for is the dependency in the flow 
direction, x. 

The biggest advantages of the ID model are the relative ease of implementa- 
tion, and comparative low computational cost in execution. Therefore, the use of 
full mult i- dimensional flow modeling, assuming its viability within the available 
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computational resources, is justified only when the ID model fails to capture the 
essential physical picture of the flow system. However, there are several limitations 
and disadvantages of the ID model that restrict its use. These limitations include, 
among other things, the Newtonian assumption, simplified fiow geometry and the 
one- dimensional nature. 

2 Theoretical Background 

The widely-used one-dimensional form of the Navier-Stokes equations to describe 
the fiow in a vessel; assuming laminar, incompressible, axi-symmetric, Newtonian, 
fully-developed fiow with negligible gravitational body forces; is given by the follow- 
ing continuity and momentum balance relations with suitable boundary conditions 
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In these equations, A is the vessel cross sectional area, t is the time, Q is the 
volumetric fiow rate, x is the axial coordinate along the vessel, L is the length of 
the vessel, a (= ■' ^--^ with u and u being the fiuid local and mean axial speed 
respectively) is the momentum fiux correction factor, p is the fiuid mass density, p 
is the local pressure, and /t is a viscosity friction coefficient which is usually given 
by K = 2TTai'/{a — 1) with u being the fiuid kinematic viscosity defined as the 
ratio of the dynamic viscosity fi to the mass density. These equations supported 
by appropriate compatibility and matching conditions are used to describe the ID 
fiow in a branched network of vessels. The equations, being two in three variables, 
Q A and p, are normally coupled with the following pressure-area relation in a 
distensible vessel to close the system and obtain a solution 
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p = Po + f{A) 



(3) 



In this relation, p and Po are the local and reference pressure respectively, and f{A) 
is a function of area which may be modeled by the following relation 
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In these equations, Ao and ho are respectively the vessel cross sectional area and wall 
thickness at reference pressure po-, while E and <;^ are the Young's elastic modulus 
and Poisson's ratio of the vessel wall. Similar variants of this ID flow model 
formulation can also be found in the literature (see for example [7, 10, 18, 19]). 

The continuity and momentum equations are usually casted in matrix form [2, 
11, 12, 16] which is more appropriate for numerical manipulation and discretization. 
In matrix form these equations are given by 
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It should be remarked that the second term of the second row of the F matrix 
can be obtained from the third term of the original momentum equation as follow 
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3 Weak Form of ID Flow Equations 
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On multiplying Equation 6 by weight functions and integrating over the solution 
domain, x, the following is obtained 
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where Q is the solution domain, and a; is a vector of arbitrary test functions. On 
integrating the second term of Equation 10 by parts, the following weak form of 
the preceding ID flow system is obtained 



/ -—- ■ udx - F ■ -^dx + / B • ujdx + [F ■ ojUq = 
in ot J^ dx in 



:iii 



where dfl is the boundary of the solution domain. This weak formulation, cou- 
pled with suitable boundary conditions, can be used as a basis for flnite element 
implementation in conjunction with an iterative scheme such as Newton-Raphson 
method. 
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4 Finite Element Solution 

There are two major cases to be considered in the finite element solution of the 
stated fiow problem: single vessel and branched network where each one of these 
cases can be time- independent or time-dependent. These four cases are outlined in 
the following three subsections. 

4.1 Single Vessel Time-Independent Flow 

The single vessel time-independent model is based on dropping the time term in the 
continuity and momentum governing equations to obtain a steady-state solution. 
This should be coupled with pertinent boundary and compatibility conditions at 
the vessel inlet and outlet. The details are given in the following. 

In discretized form, Equation 11 without the time term can be written for each 
node Ni{Ai, Qi) as 
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fi 
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where R is a vector of the weak form of the residuals and 
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where q is an index for the A''„ quadrature points, ( is the quadrature point coor- 
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dinate and 
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(15) 



with n being the number of nodes in a standard element. Because of the non- 
hnear nature of the problem, an iteration scheme, such as Newton-Raphson, can 
be utilized to construct and solve this system of equations based on the residual. 
The essence of this process is to solve the following equation iteratively and update 
the solution until a convergence criterion based on reaching a predefined error 
tolerance is satisfied 



JAU 
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In this equation, J is the jacobian matrix, AU is the perturbation vector, and 
R is the weak form of the residual vector. For a vessel with n nodes, the Jacobian 
matrix, which is of size 2n x 2n, is given by 
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where the subscripts stand for the node indices, while the vector of unknowns, 
which is of size 2n, is given by 
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In the finite element implementation, the Jacobian matrix is usually evaluated 
numerically by finite differencing, i.e. 
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AAi AQi 

The procedure to obtain a solution is summarized in the following scheme 

1. Start with initial values for Ai and Qi in the U vector. 

2. The system given by Equation 16 is constructed where the weak form of the 
residual vector R may be calculated in each iteration / (= 0, 1, . . . , M) as 



R/ 



/n(UO 
^n(UO 

3. The jacobian matrix is calculated from Equation 19. 
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(20) 



4. System 16 is solved for AU, i.e. 
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AU = -J-^R (21) 

5. U is updated to obtain a new U for the next iteration, that is 

U,+i = U, + AU (22) 

6. The norm of the residual vector is calculated from 



m = V'" + '" + --- + 'l (23) 

where e^ is the ith entry of the residual vector and A^ (= 2n) is the size of 
the residual vector. 

7. This cycle is repeated until the norm is less than a predefined error tolerance 
(e.g. 10^^) or a certain number of cycles is reached without convergence. In 
the last case, the operation will be aborted due to failure and may be resumed 
with improved finite element parameters. 

With regard to the boundary conditions (BC), two types of Dirichlet conditions 
can be applied: pressure and volumetric flow rate, that is 

A - Abc = (for area BC) & Q- Qbc = (for flow BC) (24) 

These conditions are imposed by replacing the residual function of one of the 
governing equations (the continuity equation in our model) for the boundary nodes 
with one of these constraints. 

Imposing the boundary conditions as constraints in one of the two governing 
equations is associated with imposing compatibility conditions, arising from pro- 
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jecting the differential equations in the direction of the outgoing characteristic vari- 
ables [20] , at the inlet and outlet by replacing the residual function contributed by 
the other governing equation with these conditions. The compatibility conditions 
are given by 
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where H is the matrix of partial derivative of F with respect to U, while the 
transposed left eigenvectors of H are given by 
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Hence, Equation 25 reduces to 
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that is 
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(30) 
In the last relation, the minus sign is used for the inflow boundary while the plus 
sign for the outflow boundary. The compatibility conditions, given by Equation 
30, replace the momentum residual at the boundary nodes. 

4.2 Single Vessel Time-Dependent Flow 

The aforementioned time-independent formulation can be extended to describe 
transient states by including the time terms in the residual equations in association 
with a numerical time-stepping method such as forward Euler, or backward Euler 
or central difference. The time-dependent residual will then be given (in one of the 
aforementioned schemes) by 

R^+^^* = / ^ ■ ^dx + R^+/* = (31) 

where R is the weak form of the residual, TD stands for time-dependent and TI 
for time- independent. The time-dependent jacobian follows 



-rt+At _ ^r^TD (nc)\ 
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Again, we have 



AU = -J-^R (33) 



and 



U,+i = U, + AU (34) 
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where the symbols represent time-dependent quantities and / represents Newton 
iterations. 

With regard to the boundary nodes, a steady-state or time-dependent bound- 
ary conditions could be applied depending on the physical situation while a time- 
dependent compatibility conditions should be employed by adding a time term to 
the time-independent compatibility condition, that is 



Ctd — h,2~^ + ^TI — 



(35) 



where Cti is the time-independent compatibility condition as given by Equation 
30, while Ctd is the time-dependent compatibility condition, that is 
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where the time derivatives can be evaluated by finite difference, e.g. 
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A sign convention similar to that outlined previously should be followed. An algo- 
rithmic code of the time-dependent module is presented in Algorithm 1. 



4.3 Branched Network 

To extend the time-independent and time-dependent single vessel model to time- 
independent and time-dependent branched network of interconnected vessels, match- 
ing constraints at the branching nodes are required. These nodes are treated as 
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Initialize time: t = t^ 

Initialize U*° 

for j ^— 1 to numberO fTimeSteps do 

Increment time by At 

for i ^ 1 to MaximumNumberO f Newtonlterations do 
Find R^+,^* = 4 y!i^ . ojdx + R^+/* = 

Find J*+^* - ^^'^'''^^^ 

Find AU*+^* = - (J-iR)™ * 
Find U*+f * = U*+^* + AU*+^* 
Update: U*+^* = U*+f * 
if (convergence condition met) then 

Exit loop 
else 

if {MaximumNumberO f Newtonlterations reached) then 
Declare failure 
Exit 
end if 
end if 
end for 

Update: U* = U*+^* 
end for 
Solution: U*+^* 
End 



Algorithm 1: Algorithmic code for the time-dependent module. 

discontinuous joints where each segment connected to that junction has its own 
index for that junction although they are spatially identical. The matching con- 
straints are derived from the conservation of flow rate for incompressible fluid, and 
the Bernoulli energy conservation principle for inviscid flow. More specifically, at 
each n-segment branching node, n distinctive constraints are imposed: one repre- 
sents the conservation of flow which involves all the segments at that junction, while 
the other (n — 1) constraints represent the Bernoulli principle with each Bernoulli 
constraint involving two distinctive segments. These constrains are summarized in 
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the following relations 



and 



Y,Qi = ^ (39) 



i=l 



Pk + ^pul -pi- -pu^ = (40) 

where k and / are indices of two distinct segments, and u (= ^) is the fluid speed 
averaged over the vessel cross section. In Equation 39 a directional flow is assumed 
by attaching opposite signs to the inflow and outflow. The matching constraints, 
which replace the residuals of one of the governing equations (continuity), are cou- 
pled with compatibility conditions, similar to the ones used for the single vessel, 
where these conditions replace the residual of the other governing equation (mo- 
mentum). The sign convention for these compatibility conditions should follow 
the same rules as for the boundary conditions, that is minus sign for inflow and 
plus sign for outflow. This branching model can be applied to any branching node 
with connectivity n > 2. The special case of n = 2 enables flexible modeling 
of discontinuous transition between two neighboring segments with different cross 
sectional areas. Suitable pressure or flux boundary conditions (which for the time- 
dependent case could be time-independent, or time-dependent over the whole or 
part of the time stepping process) should also be imposed on all boundary nodes 
of the network. With regard to the other aspects of the time-independence and 
time-dependence treatment, the network model should follow the same rules as for 
single vessel time-independent and time-dependent models which are outlined in 
the previous sections. 
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5 Non-Dimensionalized Form 

To improve convergence, the aforementioned dimensional forms of the governing, 
boundary, compatibihty, and matching equations can be non-dimensionahzed by 
carefully-chosen scale factors. The following scale factors are commonly used to 
scale the model parameters: 

Q ~ nRlUo A^TiRl p -~^ pU^ x ~ A ^ ~ 7^ (41) 

where Ro, Uo, and A are respectively typical values of the radius, velocity and length 
for the flow system. In the following we demonstrate non-dimensionalization of the 
flow equations by a few examples followed by stating the non-dimensionalized form 
of the others. 

5.1 Non-Dimensionalized Navier-Stokes Equations 

Continuity equation 1st form (Equation 1): 



dA dQ 



aj^^aj^r^Wl^^ (43) 

that is 



dt' dx' 
where the prime indicates a non-dimensionalized value 

Continuity equation 2nd form (Equation 6): 



^■^''^«' (44) 
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Same as Equation 44. 
Momentum equation 1st form (Equation 2): 
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dt' '^dx'\A'J ' dx' ' {a-l)RlUoA' 
Momentum equation 2nd form (Equation 6): 



dQ' d (Q"\ , A'dp' , 2avX Q' 



5.2 Non-Diniensionalized Compatibility Condition 20 



f-im-3^1---!- 






A(9t' \7cRldx' \ A' J A3pA'„ 9a;' 7ri?2A' 



\dt' Xdx' \ A' J X3pA'^ dx' A' 



UoQ' ^ jU^Q''^ ,„2 ^A , ^' l^ \ ^°^Q' 






Time-dependent term of compatibility condition: 
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5.2 Non-Dimensionalized Compatibility Condition 

Time-independent term of compatibility condition: 
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5.3 Non-Dimensionalized Matching Conditions 

Flow conservation: 

Q'l - Q'2 - Qs = (58) 

Bernoulli: 



P'k + 2^'fc ~P'i^ 2^'^i = (^9) 



5.4 Non-Dimensionalized Boundary Conditions 

A' - A'bc = (for area BC) & Q' - Q'bc = (for Sow BC) (60) 

6 Validation 

The different modules of the the ID finite element flow model are validated as 
follow 

• Time-independent single vessel: the numeric solution should match the ana- 
lytic solution as given by Equation 63 which is derived in Appendix A. Also, 
the boundary conditions should be strictly satisfied. 

• Time-dependent single vessel: the solution should asymptotically converge 
to the analytic solution on imposing time-independent boundary conditions. 
Also, the boundary conditions should be strictly satisfied at all time steps. 

• Time-independent network: four tests are used to validate the numeric so- 
lution. First, the boundary conditions should be strictly satisfied. Second, 
the conservation of mass (or conservation of volume for incompressible flow), 
as given by Equation 39, should be satisfied at all branching nodes (bridge. 
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bifurcation, trifurcation, etc.). A consequence of this condition is that the 
sum of the boundary inflow (sum of Q at inlet boundaries) should be equal 
to the sum of the boundary outflow (sum of Q at outlet boundaries). Third, 
the conservation of energy (Bernoulli's principle for inviscid flow), as given by 
Equation 40, should be satisfied at all branching nodes. Fourth, the analytic 
solution for time-independent flow in a single vessel, as given by Equation 63 
in Appendix A, should be satisfied by all vessels in the network with possible 
exception of very few vessels with odd features (e.g. those with distorted 
shape such as extreme radius-to-length ratio, and hence are susceptible to 
large numerical errors). The fourth test is based on the fact that the sin- 
gle vessel solution is dependent on the boundary conditions and not on the 
mechanism by which these conditions are imposed. 

• Time-dependent network: the solution is validated by asymptotic conver- 
gence to the time- independent solution, as validated by the four tests outlined 
in the previous item, on imposing time-independent boundary conditions. 

The solutions may also be tested qualitatively by static and dynamic visual- 
ization for time-independent and time-dependent cases respectively to verify their 
physical sensibility. Other qualitative tests, such as comparing the solutions of 
different cases with common features, may also be used for validation. 

It should be remarked that Equation 63 contains three variables: x A and Q, 
and hence it can be solved for one of these variables given the other two. Solving 
for A and Q requires employing a numeric solver, based for example on a bisection 
method; hence the best option is to solve for x and compare to the numeric solution. 
This in essence is an exchange of the role of independent and dependent variables 
which has no effect on validation. Alternatively, Equation 76 can be used to verify 
the solution directly by using the vessel inlet and outlet areas. In fact Equation 76 
can be used to verify the solution at any point on the vessel axis by labeling the 
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area at that point as Aou, as explained in Appendix A. 

7 General Notes 

7.1 Implementation 

• The model described in this report was implemented and tested on both sin- 
gle vessels and networks of vessels for time-independent and time-dependent 
cases and it produced valid results. The implementation is based on a 
Galerkin method, where the test functions are obtained from the same space 
as the trial basis functions used to represent the state variables, with a La- 
grange interpolation associated with a Gauss quadrature integration scheme 
(refer to Appendix B for Gauss quadrature tables). Many tests have been 
carried out to verify various aspects of the ID model. These tests involved 
many synthetic and biological networks which vary in their size, connectivity, 
number and type of branching nodes, type of meshing, and so on. The tests 
also included networks with and without loops although the great majority 
of the networks contain loops. Some of the networks involved in these tests 
consist of very large number of vessels in the order of hundreds of thousands 
with much more degrees of freedom. Non-dimensional form, as well as dimen- 
sional form, was tested on single vessels and branched networks; the results, 
after re-scaling, were verified to be identical to those obtained from the di- 
mensional form. The checks also included h and p convergence tests which 
demonstrated correct convergence behavior. 

• To be on the safe side, the order of the quadrature should be based on the sum 
of orders of the interpolating functions, their derivatives and test functions. 
The adopted quadrature order scheme takes the highest order required by 
the terms. 



7.2 Solution 24 

• A constant delta may be used in the evaluation of the Jacobian matrix by 
finite difference. A suitable value for delta may be AAi = AQi = 10~^ or 
10-8. 

7.2 Solution 

• Negative flow in the solution means the flow direction is opposite to the vessel 
direction as indicated by the vessel topology, that is the flow rate of a segment 
indexed as A^^i N2 will be positive if the flow is from A^^i to N2 and negative 
if the flow is from N2 to Ni. 

• Interpolation polynomials of various degrees (p) in association with differ- 
ent meshing (h) should be used to validate the convergence behavior. The 
convergence to the correct solution should improve by increasing p and de- 
creasing h. L^ error norm may be used as a measure for convergence; it is 
given by 

\ 1/2 

{Sa-Snfdx] (61) 

X J 

where Sa and Sn are the analytic and numeric solutions respectively, and X is 
the solution domain. The integration can be performed numerically using, for 
example, trapezium or Simpson's rules. The error norm should fall steadily 
as h decreases and p increases. 

With regard to the previously outlined implementation of the ID model (see 
§ 7.1), typical solution time on a typical platform (normal laptop or desktop) 
for a single time-independent simulation on a typical ID network consisting 
of hundreds of thousands of degrees of freedom is a few minutes. The flnal 
convergence is normally reached within 5-7 Newton iterations. The solution 
time of a single time step for the time-dependent case is normally less than 
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the solution time of the equivalent time-independent case, and the number of 
the required Newton iterations of each time step in the time-dependent case 
is normally less than that for the corresponding time-independent case. 

• Since there are many sources of error and wrong convergence, each acquired 
solution should be verified by the aforementioned validation tests (see § 6). 
The ID finite element code should be treated as a device that suggests solu- 
tions which can be accepted only if they meet the validation criteria. 

7.3 Non-Dimensionalization 

• On implementing the non-dimensionalized form (as given in § 5) in the finite 
element code, all the user needs is to scale the primed input values either 
inside or outside the code; the results then should be scaled back to obtain 
the dimensionalized solution. 

• Different length scales can be utilized as long as they are in different orien- 
tations (e.g. vessel length and vessel radius) and hence linearly independent; 
otherwise the physical space will be distorted in non-systematic way and 
hence may not be possible to restore by scaling back. 

7.4 Convergence 

A number of measures, outlined in the following points, can speedup convergence 
and help avoiding convergence failure. 

• Non-dimensionalization which requires implementation in the finite element 
code (as given in § 5) where the input data is non-dimensionalized and the 
results are re-dimensionalized back to the physical space. 

• Using different unit systems, such as m.kg.s or mm.g.s or m.g.s, for the input 
data and parameters. 
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• Scaling the model up or down to obtain a similarity solution which can then 
be scaled back to obtain the final results. 

• Increasing the error tolerance of the solver for convergence criterion. How- 
ever, the use of relatively large error tolerance can cause wrong convergence 
and hence should be avoided. It may be recommended that the maximum 
allowed error tolerance for obtaining a reliable solution must not exceed 10~^. 
Anyway, the solution in all cases should be verified by the validation tests 
(see § 6) and hence it must be rejected if the errors exceed acceptable limits. 

• For time-dependent cases, the required boundary condition value can be im- 
posed gradually by increasing the inlet pressure, for instance, over a number 
of time steps to reach the final steady state value. 

• The use of smaller time steps in the time-dependent cases may also help to 
avoid convergence failure. 

It should be remarked that the first three strategies are based on the same 
principle, that is adjusting the size of the problem numbers to help the solver to 
converge more easily to the solution. 

7.5 Boundary Conditions 

• Dirichlet type boundary conditions are usually used for imposing fiow rate 
and pressure boundary conditions. The previous formulation is based on this 
assumption. 

• Pressure boundary conditions are imposed by adjusting the inlet or outlet 
area where p and A are correlated through Equation 3. 

• While pressure boundary conditions can be imposed on both inlet and outlet 
boundaries simultaneously, as well as mixed boundary conditions (i.e. inlet 
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pressure with outlet flow or inlet flow with outlet pressure or mixed on one 
or both boundaries), it is not possible to impose flow boundary conditions 
on all inlet and outlet boundaries simultaneously because this is either a 
trivial condition repeating the condition of flow conservation (i.e. Equation 
39) at the branching junctions if the inflow is equal to the outflow or it is a 
contradiction to the flow conservation condition if the inflow and outflow are 
different, and hence no solution can be found due to ambiguity and lack of 
constraints in the first case and to inconsistency in the second case. 

• Zero Q boundary condition can be used to block certain inlet or outlet vessels 
in a network for the purpose of emulating a physical situation or improving 
convergence when the blockage does not affect the solution significantly. 

• In some biological flow conditions there are no sufficient data to impose real- 
istic pressure boundary conditions that ensure biologically sensible flow in the 
correct direction over the whole vascular network. In such situations a back 
flow may occur in some branches which is physically correct but biologically 
incorrect. To avoid this situation, an inlet pressure boundary condition with 
outlet flow boundary conditions where the total outflow is split according to 
a certain physical or biological model (such as being proportional to the area 
squared) can be used to ensure sensible flow in the right direction over the 
whole network. The total amount of the outflow can be estimated from the 
inlet flow which is usually easier to estimate as it normally comes from a sin- 
gle (or few) large vessel. This trick may also be applicable in some physical 
circumstances. 

• Use may be made of an artificial single inlet boundary to avoid lack of knowl- 
edge about the pressure distribution in a multi-inlet network to ensure correct 
flow in the right direction. The inlets can be connected to a single artificial 
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node (e.g. located at their centroid) where the radii of the connecting ar- 
tificial vessels is chosen according to a physical or biological model such as 
Murray's law. This node can then be connected through a single artificial 
vessel whose radius can be computed from a physical or biological model and 
whose length can be determined from a typical L/R ratio such as 10. The 
inlet of this vessel can then be used to impose a single p ot Q biologically- 
sensible boundary condition. It should be remarked that Murray's biological 
law is given by 



m 



H'l (62) 



where r^ is the radius of the mother vessel, r^^. is the radius of the iih. daughter 
vessel, n is the number of daughter vessels which in most cases is 2, and 7 
is a constant index which according to Murray is 3, but other values like 2.1 
and 2.2 are also used in the literature. 

• Time-dependent boundary conditions can be modeled by empirical signals 
(e.g. obtained from experimental data) or by closed analytical forms such as 
sinusoidal. 

7.6 Initial Conditions 

• The convergence usually depends on the initial values of area and flow rate. 
A good option for these values is to use unstressed area with zero flow for 
start. 

7.7 Miscellaneous 

• Apart from the interpolation nodes, there are two main types of nodes in 
the finite element network: segment nodes and finite element discretization 
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nodes. The connectivity of the second type is always 2 as these nodes connect 
two elements; whereas the connectivity of the first can be 1 for the bound- 
ary nodes, 2 for the bridge nodes connecting two segments, or > 3 for the 
branching nodes (bifurcation, trifurcation, etc.). The mass and energy con- 
servation conditions can be extended to include all the segment nodes with 
connectivity > 1 by including the bridge nodes. 

• For networks, the vessel wall thickness at reference pressure, ho, can be a 
constant or vary from vessel to vessel depending on the physical or biological 
situation. Using variable thickness is more sound in biological context where 
the thickness can be estimated as a fraction of the lumen or vessel radius. 
A fractional thickness of 10-15% of the radius is commonly used for blood 
vessels [6, 16, 20-26]. For more details, refer to Appendix D. 

• The previous finite element formulation of the ID model for single vessels and 
networks works for constant-radius vessels (i.e. with constant Ao) only and 
hence to extend the formulation to variable-radius vessels the previous matrix 
structure should be reshaped to include the effect of tapering or expanding 
of the vessels. However, the vessels can be straight or curved. The size of the 
vessels in a network can also vary significantly from one vessel to another as 
long as the ID flow model assumptions (e.g. size, shape, etc.) do apply on 
each vessel. 

• The networks used in the ID flow model should be totally connected, that 
is any node in the network can be reached from any other node by moving 
entirely inside the network vessels. 

• Different time stepping schemes, such as forward or backward Euler or central 
difference, can be used for implementing the time term of the time-dependent 
single vessel and network modules although the speed of convergence and 
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quality of solution vary between these schemes. The size of the time step 
should be chosen properly for each scheme to obtain equivalent results from 
these different schemes. 

• Although the ID model works on highly non-homogeneous networks in terms 
of vessels length without discretization, a scheme of homogeneous discretiza- 
tion may be employed by using a constant element length, h, over the whole 
network as an approximation to the length of the discretized elements. The 
length of the elements of each vessel is then obtained by dividing the ves- 
sel evenly to an integer number of elements with closest size to the given 
h. Although discretization is not a requirement, since the ID model works 
even on non-discretized networks, it usually improves the solution. More- 
over, discretization is required for obtaining a detailed picture of the pressure 
and flow fields at the interior points. Use of interpolation schemes higher 
than linear (with and without discretization) also helps in refining the vari- 
able fields. Also, for single vessel the solution can be obtained with and 
without discretization; in the first case the discretized elements could be of 
equal or varying length. The solution, however, should generally improve by 
discretization. 

• Although the ID model works on non-homogeneous networks in terms of 
vessels radius, an abrupt transition from one vessel to its neighbor may hinder 
convergence. 

• In general, the time-dependent problem converges more easily than its equiva- 
lent time-independent problem. This may be exploited to obtain approximate 
time-independent solutions in some circumstances from the time- dependent 
module as the latter asymptotically approaches the time-independent solu- 
tion. 



7. 7 Miscellaneous 3 1 

• The correctness of the solutions mathematically may not guarantee physio- 
logical, and even physical, sensibility since the network features, boundary 
conditions, and model parameters which in general highly affect the flow pat- 
tern, may not be found normally in real biological and physical systems. The 
quality of any solution, assuming its correctness in mathematical terms, de- 
pends on the quality of the underlying model and how it reflects the physical 
reality. 

• Because the ID model depends on the length of the vessels but not their 
location or orientation, a ID coordinate system, as well as 2D or 3D, can 
be used for coordinating the space. The vessels can be randomly oriented 
without effecting the solution. A multi-dimensional space may be required, 
however, for consistent and physically-correct description of the networks. 

• The reference pressure, po, in Equation 3 is usually assumed zero to simplify 
the relation. 
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8 Conclusions 

The one-dimensional Navier-Stokes formulation is widely used as a realistic model 
for the flow of Newtonian fluids in large vessels with certain simplifying assump- 
tions, such as axi-symmetry. The model may also be coupled with a pressure-area 
constitutive relation and hence be extended to the flow in distensible vessels. Nu- 
merical implementation of this model based on a finite element method with suit- 
able boundary conditions is also used to solve the time-independent and transient 
flow in single vessels and networks of interconnected vessels where in the second 
case compatibility and matching conditions, which include conservation of mass 
and energy, at branching nodes are imposed. Despite its comparative simplicity, 
the ID flow model can provide reliable solutions, with relatively low computational 
cost, to many flow problems within its domain of validity. The current document 
outlined the analytical and numerical aspects of this model with theoretical and 
technical details related to implementation, performance, methods of improvement, 
validation, and so on. 
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Nomenclature 

a correction factor for axial momentum flux 

P parameter in the pressure-area relation 

7 Murray's law index 

e residual error 

( quadrature point coordinate 

K viscosity friction coefficient 

H fluid dynamic viscosity 

u fluid kinematic viscosity {u = -) 

p fluid mass density 

^ Poisson's ratio of vessel wall 

ip basis function for flnite element discretization 

LO vector of test functions in the weak form of flnite element 

Q solution domain 

dQ boundary of the solution domain 

A vessel cross sectional area 

Abc boundary condition for vessel cross sectional area 

Ain vessel cross sectional area at inlet 

Ao vessel cross sectional area at reference pressure 

B matrix of force terms in the ID Navier-Stokes equations 

E Young's modulus of vessel wall 

f{A) function in pressure-area relation 

F matrix of flux quantities in the ID Navier-Stokes equations 
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h length of element 

ho vessel wall thickness at reference pressure 

H matrix of partial derivative of F with respect to U 

J Jacobian matrix 

L length of vessel 

9T norm of residual vector 

p local pressure 

p order of interpolating polynomial 

Po reference pressure 

q dummy index for quadrature point 

Q volumetric flow rate 

Qbc boundary condition for volumetric flow rate 

r radius 

R weak form of residual vector 

Sa analytic solution 

Sn numeric solution 

t time 

At time step 

u local axial speed of fluid 

u mean axial speed of fluid 

U vector of finite element variables 

AU vector of change in U 

X vessel axial coordinate 
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9 Appendix A: Derivation of Time-Independent 
Analytical Solution for Single Vessel 

The following analytical relation linking vessel axial coordinate x to cross sectional 
area A, cross sectional area at inlet Ain, and volumetric flow rate Q for time- 
independent flow can be derived and used to verify the finite element solution 

aQHniA/A,^)-^(A^/^-Ar) 
X = ^ L (63) 

The derivation is outlined in the following. For time-independent flow, the 
system given by Equation 6 in matrix form, will become 



?^ = X G [0, /] , t > (64) 

ox 



TM^^f")^^^-^ -l»''l-*^° 



dx 
that is Q as a function of x is constant and 



(65) 



a (-Q\^M?A^S=, (66) 



dA\ A 2,pAo J dx A 



I.e. 



"«^ + -l^M ^ + 4 = (67) 



A^ 2pAo J dx A 
which by algebraic manipulation can be transformed to 

dA ~ -kQ 

On integrating the last equation we obtain 



(6^ 
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X = — -^ + C (69) 

where C is the constant of integration which can be determined from the boundary 
condition at x = with A = Ain, that is 

C = ^ ^^^° 70 

i.e. 



aQ2ln(AMi„)-JL(^5/2_4/2^ 



^= #^ '- P^' 

For practical reasons, this relation can be re-shaped and simplify to reduce the 
number of variables by the use of the second boundary condition at the outlet, as 
outlined in the following. When x = L, A = Aou where L is the vessel length and 
Aou is the cross sectional area at the outlet, that is 



aQ^ In (A,„M,„) - 5^ I^AlL' - A^ 
which is a quadratic polynomial in Q i.e. 



' ^Q-^ '- (^^) 



a In (^„M,„) Q' + nLQ + ^^ (a^/^ - A^i') = (73) 



a In {A^JA,^) Q' + nLQ + ^^ I^A^' - A^i') = (74) 

with a solution given by 



-kL ± Jk^L^ - 4a In (A.JAou) 5^ [aH' - A^^ 
^ ^ 2aln(A,„M,„) ^^^^ 
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which is necessarily real for Ain > Aou which can always be satisfied for normal 
fiow conditions by proper labeling. For a fiow physically-consistent in direction 
with the pressure gradient, the root with the plus sign should be chosen, i.e. 



-kL + Jk^L^ - 4a In (A^JAou) ^^ [aIL' - A^ 
^ ^ 2«ln(A,„M„J ^^^^ 

This, in essence, is a relation between fiow rate and pressure drop (similar to 
the Hagen-Poiseuille law for rigid vessels) although for elastic vessels the fiow rate, 
as given by Equation 76, does not depend on the pressure difference (as for rigid 
vessels) but on the actual inlet and outlet pressures. 

Although Equation 76 may look a special case of Equation 63 as it involves only 
the vessel two end areas, Aou may be assumed to be the area at any point along 
the vessel axis, with L being the distance form the vessel inlet to that point, and 
hence this relation can be used to verify the finite element solution at any point on 
the vessel. 
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10 Appendix B: Gauss Quadrature 

In this appendix we list points and weights of Gauss quadrature for polynomials 
of order 1-10 which may not be easy to find. 



I— 1 





Table 1: Gauss quadrature points and weights for polynomial order, p, 1-10 assuming a 0-1 master element. 




p 


Points 


1 

2 
3 

4 
5 
6 
7 
8 
9 
10 


0.500000000000000 
0.211324865405188 
0.112701665379259 
0.069431844202974 
0.046910077030669 
0.033765242898424 
0.025446043828621 
0.019855071751232 
0.015919880246187 
0.013046735741414 


0.788675134594812 
0.500000000000000 
0.330009478207572 
0.230765344947159 
0.169395306766867 
0.129234407200303 
0.101666761293187 
0.081984446336682 
0.067468316655508 


0.887298334620741 
0.669990521792428 
0.500000000000000 
0.380690406958402 
0.297077424311301 
0.237233795041836 
0.193314283649705 
0.160295215850488 


0.930568155797026 
0.769234655052841 
0.619309593041598 
0.500000000000000 
0.408282678752175 
0.337873288298095 
0.283302302935377 


0.953089922969331 
0.830604693233133 
0.702922575688699 
0.591717321247825 
0.500000000000000 
0.425562830509185 


0.966234757101576 
0.870765592799697 
0.762766204958164 
0.662126711701905 
0.574437169490815 


0.974553956171380 
0.898333238706813 
0.806685716350295 
0.716697697064623 


0.980144928248768 
0.918015553663318 
0.839704784149512 


0.984080119753813 
0.932531683344493 


0.986953264258586 




Weights 


1 
2 
3 

4 
5 
6 
7 
8 
9 
10 


1.000000000000000 
0.500000000000000 
0.277777777777778 
0.173927422568727 
0.118463442528095 
0.085662246189585 
0.064742483084435 
0.050614268145188 
0.040637194180787 
0.033335672154344 


0.500000000000000 
0.444444444444444 
0.326072577431273 
0.239314335249683 
0.180380786524069 
0.139852695744638 
0.111190517226687 
0.090324080347429 
0.074725674575291 


0.277777777777778 
0.326072577431273 
.284444444444444 
0.233956967286345 
0.190915025252559 
0.156853322938943 
0.130305348201467 
0.109543181257991 


0.173927422568727 
0.239314335249683 
0.233956967286345 
0.208979591836735 
0.181341891689181 
0.156173538520001 
0.134633359654998 


0.118463442528095 
0.180380786524069 
0.190915025252559 
0.181341891689181 
0.165119677500630 
0.147762112357376 


0.085662246189585 
0.139852695744638 
0.156853322938943 
0.156173538520001 
0.147762112357376 


0.064742483084435 
0.111190517226687 
0.130305348201467 
0.134633359654998 


0.050614268145188 
0.090324080347429 
0.109543181257991 


0.040637194180787 
0.074725674575291 


0.033335672154344 



to 

Co 

s 

a 
§ 
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11 Appendix C: Biological Parameters 

In this appendix, we suggest some biologically-realistic values for the ID flow model 
parameters in the context of simulating blood flow in large vessels. 

1. Blood mass density (p): 1050 kg.m'^ [6, 7, 16, 25, 27-31]. 

2. Blood dynamic viscosity (/i): 0.0035 Pa.s [4, 6, 7, 15, 16, 20, 25, 27, 30, 32-34]. 

3. Young's elastic modulus (E): 100 kPa [4, 6, 13, 15, 16, 20, 22, 25, 27, 34-38]. 
Also see [25, 39] on shear modulus. 

4. Vessel wall thickness [ho): this, preferably, is vessel dependent, i.e. a fraction 
of the lumen or vessel radius according to some experimentally-established 
mathematical relation. The relation between wall thickness and vessel inner 
radius is somehow complex and vary depending on the type of vessel (e.g. 
artery or capillary). For arteries, the typical ratio of wall thickness to inner 
radius is about 0.1-0.15, and this ratio seems to go down in the capillaries and 
arterioles. Therefore a typical value of 0.1 seems reasonable [6, 16, 20-26]. 

5. Momentum correction factor (a): 4/3 = 1.33 assuming Newtonian flow. A 
smaller value, e.g. 1.2, may be used to account for non-Newtonian shear- 
thinning effects [2, 3, 7, 16, 17, 20, 28, 40]. 

6. Time step (At): 1.0-0.1 ms [4, 7, 9, 10, 13, 15, 16, 18-20, 25, 29, 34, 41-45]. 

7. Pressure step (Ap): 1.0-5.0 kPa [4, 7, 16, 29]. 

8. Time of heart beat: 0.85 s assuming 70 beats per minute. 

9. Poisson ratio (^): 0.45 [2, 4, 6, 13, 22, 25, 27, 34, 37, 39, 41]. 



